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TRIANGULAR  FACTORIZATION  AND  INVERSION 
BY  FAST  MATRIX  MULTIPLICATION 

James  R.  Bunch 
and 

John  E.  Hopcroft 


1.  INTRODUCTION 

Strassen  [3]  has  given  an  algorithm  using  non-commutative 

multiplication  which  computes  the  product  of  two  matrices  of 

order  2  by  7  multiplications  and  18  additions.  Thus  the  product 

k  3  k 

of  two  matrices  of  order  m2  could  be  computed  by  m  7  multi- 

2  k  k  2 

plications  and  (5+m)m  7  -  6  (m2  )  additions. 

Strassen  uses  block  LDU  factorization  (Householder  [2] , 

p.  126)  recursively  to  compute  the  inverse  of  a  matrix  of  order 

m2k  by  m2k  divisions,  £  |m37k  -  m2k  multiplications,  and 
6  2  k  k  2 

<  |-(5+m)m  7  -  7  (m2  )  additions.  The  inverse  of  a  matrix  of 

order  n  could  then  be  computed  by  £  5.64nAo^2^  arithmetic 
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if  A^  and  h  are  non -singular. 

Since  the  algorithm  is  applied  recursively,  it  will  fail 
whenever  the  inversion  of  a  singular  principal  submatrix  in  any 
of  the  reduced  matrices  is  required. 

For  example,  the  block  LDU  factorization  fails  to  exist 
for  a  matrix  as  simple  as 


0  0  0  1  | 

0  0  10 

0  10  0 

10  0  0 
»  _ 

Every  principal  submatrix  in  every  reduced  matrix  is  non¬ 
singular  if  A  is  symmetric  positive  definite,  strictly  diagonally 
dominant,  or  irreducibly  diagonally  dominant  (Varga  [4],  p.  23). 
However,  if  A  is  only  non-singular,  then  we  must,  in  general, 
pivot  (i.e.,  interchange  rows  or  columns)  in  order  to  obtain  a 
(point  or  block)  LDU  factorization.  If  A  is  non-singular,  then 
there  exist  permutation  matrices  P^,  Q^,  Q ^  such  that 

AP1'  ^1A'  ^2AP2  ^ave  (P°int  or  block)  LDU  factorizations  (cf. 

Forsythe  and  Moler  [1],  p.  36). 

In  Sections  2  and  3  we  show  that  by  employing  pivoting 
we  can  use  Strassen's  fast  matrix  multiplication  algorithm  to 
obtain  the  LU,  or  LDU,  decomposition  of  any  non-singular  matrix 

/ 


! 

« 

-  3  - 

Jc  £qq  1 

of  order  n  =  2  in  <  (3.64)n92  operations,  and  hence  its 
inverse  in  <  (10.18)n*og2^  operation,  where  an  operation  is 
defined  to  be  a  multiplication,  division,  addition,  or  subtrac¬ 
tion. 

In  Section  4  we  modify  the  algorithm  that  we  can  find 
triangular  factorizations  in  <  (2.04)nAo^27  operations  and 
inverses  in  <  (5.70)n^°^27  operations  when  n  =  2^.  Then, 
for  arbitrary  n,  we  can  find  triangular  factorizations  in 
<  (3.07)n*°^27  operations  and  inverses  in  <  (7.46)n^og27 
operations . 
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2.  THE  BASIC  ALGORITHM 

Jc 

For  simplicity,  let  M  be  of  order  n  =  2  with  detM^O. 

Let  =  M.  We  shall  construct  a  sequence  P1,  P2,  . ..,  Pn  1 

of  permutation  matrices  so  that  M  =  LUP,  i.e.  MP  1  =  LU,  where 

P  E  P1  P2  ***  Pn_1  is  a  permutation  matrix,  L  =  L1  L2  Ln  1 

is  unit  lower  triangular,  U  is  upper  triangular,  and 

n  i-15 

det  M  =  (det  P)  det  U  =  ±  n  u.  .  .  Since  (PA)  =  P  here. 


p-l.  o  pn-l  ^  p2  P1/  and 

M_1  =  p"1!)"3!,"1  =  pH"1.  .  .P2P1U_1(Ln"1)“1  ...  (L2)"1(L1)"1, 

where  (L1)  1  =  21  -  L^. 

We  define  the  algorithm  sequentially  for  1  £  i  <_  n-1 
as  follows. 


Let  =  { j  :  i.=l,  i=i^_^.2k  1+^jc_2^k  2+. .  ,+i^21+iQ2®} ; 


let  t  =  meix{j:j  t  Bi>  ,  s  =  min{j:j  B^},  and  r 

j  j 


ft  if  sj*t 
=  \  t-1  if  s=t  * 


Then  M 


■* 

«.i—  1 
M11 

m,i  — 1 
M12 

0 

i-1 

Mi_1 

wM21 

M22 

m 

,  where  is  a  non-singular 


upper  triangular  matrix  of  order  i-1,  ls  (i“i)x(n”i+l)  • 

0  is  the  (2r+1-i+l) x  (i-1)  zero  matrix,  M^1  is  (n-2r+1)x  (i-1) , 
M22  is  (n-i+1) x (n-i+1) ,  and  M1  is  non-singular. 


Since  2r+*-i+l  >  0  and  M1-^  is  non-singular,  there 


exists  a  non-zero  element  in  the  first  row  of  M, 


Hence 


•  •  *  1  * 

there  exists  a  permutation  matrix  P1  such  that  N1  =  M1-1?1, 
nf .  0  ,  and  N*  can  be  partitioned  as : 


N1  - 


i  e1! 


i  — :- 

°  I  G1  L  H1 

r - 


x1  I 


,  where  U*  is  (i-2s)x(i-2s)  ,  V*  is 


(i-2s)x (n-i+2s) ,E*  and  are  2sx2s,  and  H^"  are  2sx(n-i), 

0  is  the  (2r+1-i+2s)x (i— 2s )  zero  matrix,  X*  is  (n-2r+1)x ( i— 2s )  , 

i  s  s  i  i 

and  Y  is  (n-i-2  )x(n-i+2  ).  Further,  U  and  E  are  non¬ 


singular  upper  triangular. 


Let  Z1  =  G1 (E1) -1  and  L1  = 


i-2s 


O  Vi 


where  1^  is  the  identity  matrix  of  order  j. 

*  i  ■  • 

Define  M1  =  (L1)"1  N1.  Then 


M1  « 


U1  V1 

- j  T . 

I  i«  i 

o  is-j-e:. 

l2Ld 


,  where  J1  =  H1  -  Z1  F1  . 


At  the  last  step  U  =  Mn_1  is  non-singular  and  upper  tri¬ 


angular. 
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3.  OPERATION  COUNT. 


Finding  the  permutation  P*  requires  at  most  n-i  compari¬ 
sons,  and  if  PU  =  0  then  the  permutation  involves  n  element 
interchanges.  Hence  at  most  n(n-l)/2  comparisons  and  at  most 
n(n-l)  element  interchanges  are  required  to  obtain  M  =  LUp. 

The  computation  of  m”*  would  require  at  most  an  additional 
n(n-l)  element  interchanges. 

Let  an  operation  be  a  multiplication,  division,  addition, 
or  subtraction.  Let  M(n),  MT(n),  and  IT(n)  be  the  number 
of  operations  required  to  multiply  two  nxn  matrices,  to  multiply 
an  nxn  matrix  by  an  upper  triangular  nxn  matrix,  and  to  invert 

an  nxn  non-singular  upper  triangular  matrix  (we  shall  ignore 

k  k+l 

lower  order  terms).  Then  M(l)  =  1  and  M(2  )  =  7  for 

k>l.  Since  M_(2k)  =  4M_,(2k_1)  +  2M(2k~1)  +  22k-1  and 
“  11  k-1 

IT(2k)  =  2  IT(2k’1)  +  2  MT(2k"1)  ,  MT(2k)  =  2  I  4jM(2k“j"1) 


<  (ii)7 


and 


k  k_  ^  -4  k--4_l 

I_(2K)  =2  I  2J  Mt(2K  3  l 
1  j=0 


)  <  (g)7k. 


Inverting  all  the 


for 


1  <  i  <  n-1  requires 


kr-1 

>k-l  V  * 
2^ 


l 

j=0 


ItCZ3) 


<  2 


k-1 


.28. 


k-1 

V 

L 

j=0 


7vj 


<7> 


14.  „k  ,2.7 


<  (g)2 


V? 


/  \  7k 

(75)7 


operations.  Forming  all  the  multipliers  Z1  for  1  <  i  <  n-1 

k-1 

requires  2k_1  £  MT(2j)  <  (j|)  7k  operations.  Forming  all 

j=0  2 

the  reduced  matrices  J*  for  l£i<n-l  requires 


7 


k-1 

E 

3*o 


7  ,2k 

t 


(2£+l) 1 


M(23)  <  22(k-l) 


<  1  22k  (J'k  ,4'  -  ,7'  -k 


(J)  <y>  =  <T>  7  • 


k-l 

y  m(2J) 

Za  ?T~ 

3=0 


iv 

Hence,  for  n  =  2  ,  triangular  factorization  requires 

<  (§■§■)  7k  =  3.64n£ogz7  operations. 


Inverting  U  requires  <  (^-)  7k  operations  and  L-* 

k-l  .  k-l 

reauires  22k  V  M(2  }  -  (7 )  22k  V  l7)?  < 

requires  2  ^  -I^l  -  (?)  2  ^  (T>  < 

j=o 


3*o 

k3^  “  '  3; 


,7._2k  ,7.k  ,4.  ,14. _k  .  . 

(2)2  (3-)  (?)  =  (”)7  operations. 


Hence,  for  n  =  2  ,  inversion  requires  {— ^)1  <  (10. 18) n£ogz7 


operations, 


If  M  is  a  non-singular  matrix  of  order  n,  where 

2k  <  n  <  2k+^ ,  then  let  *At  -  M  0  I  .  We  can  find  the 

2  -n 

triangular  factorization  of  a  permutation  otdi,  and  hence  of  a 
permutation  of  M,  by  <  (^|-)7k+^  =  (-j^)  7k  <  (25.48)n£ogz7  , 
and  the  inverse  of  and  hence  of  M,  by  <  (-^y|)7k+^  = 

(~~)7k  <  (71.22)n£ogz7. 


4.  A  MODIFIED  ALGORITHM 


Wc  can  modify  the  algorithm  in  Section  2  so  that  the 
coefficient  of  n*°927  is  smaller  in  the  operation  counts 
of  Section  3,  In  particular,  we  find  m  and  k  such  that 
the  number  of  operations  is  minimized  subject  to  the  con- 

V 

straint  n  £  m  2  . 

First,  let  n  *  2r  »  m  2k.  Then  m  =  2r“k  =  2s  , 

and  M(2r)  <  (5+2m)m27k  =  f(s)7r  ,  where  f(s)  =  (5+2m)m27_s  = 

(5+2S+1)22s7‘s.  Since  min  f(s)  =  f(3)  =  ,  we  take 

0<s <r  * 

m  =  8,  k  =  r-3,  and  use  regular  multiplication  and  inversion 

r  197  r 

for  8x8  matrices.  Then  M(2  )  <_  p^4)7r  for  r  £  0  (rather 

than  M(2r)  £  (7)7r  in  Section  3).  Hence  each  coefficient 

1  192 

in  Section  3  is  multiplied  by 

Triangular  factorization  requires  <  <  (2 . 04) nAog2 7 

operations,  and  inversion  requires  <  (^yg-)  (jly) 7*  <  [5. 70)n^og27. 

Now  let  n  be  arbitrary.  Taking  k  =  [£og2n-4]  and 
ra  «  ln2~k]+l  (cf.  [31),  we  have  n  £  m2k  and  (5+2m)m27k  < 
(4.7)n*°927. 

k-1 

Now  MT(m2k)  <  2(5+2m)m27k”1  4)j  <  §<5+2m>m27k 

j=0 

and  IT(m2k)  <  |j-(5+2m)m27k"1  V  (|)j  <  ^(5+2m)m27k. 

49  2  k 

Triangular  factorization  thus  requires  <  yg-(5+2m)m  7  < 

(3.07)n*°^27  operations,  and  inversion  requires  <  i^-(5+2m)m27k  < 
(7.46)n£o^27  operations. 
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5.  REMARKS. 

As  seen  above,  the  coefficient  of  is  very 

sensitive  to  the  implementation  of  the  algorithm.  Another 
modification  of  the  algorithm  might  reduce  the  coefficient. 
Further,  the  bounds  we  have  given  on  the  coefficient  are 
pessimistic. 

The  algorithm  as  stated  in  Section  2  and  4  may  not 
be  numerically  stable  since  we  cannot  guarantee  that  the 
elements  in  the  reduced  matrices  are  bounded.  However, 
there  may  be  a  modification  of  our  algorithm  which  guarantees 
stability;  this  question  deserves  further  investigation. 

If  a  fast  matrix  multiplication  algorithm  were  given 
for  multiplying  two  matrices  of  order  u  in  v  multiplica¬ 
tions,  then  algorithms  similar  to  those  in  Sections  2  and  4 
could  find  the  triangular  factorization  of  a  permutation 
of  any  non-singular  matrix,  and  hence  the  inverse  of  any 
non-singular  matrix,  in  <  c  nio^uv  operations.  The  algorithms 
would  be  expressed  in  terms  of  the  expansions  of  integers 
modulo  u. 
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